library(sf)
library(dplyr)
library(ggplot2)
library(MASS) # kde2d
library(scales) # rescale
# 1) Trees in Council District 7
trees_cd7 <- trees_with_district |>
filter(council_dist == 7)
coords_cd7 <- st_coordinates(trees_cd7)
lon <- coords_cd7[, 1]
lat <- coords_cd7[, 2]
# 2) KDE on lon/lat
dens <- MASS::kde2d(
x = lon,
y = lat,
n = 150 # smaller = faster; increase later if you want smoother
)
# Turn KDE grid into a data.frame of points
grid_df <- expand.grid(
lon = dens$x,
lat = dens$y
)
grid_df$density <- as.vector(dens$z)
# Convert to sf points
density_pts <- st_as_sf(
grid_df,
coords = c("lon", "lat"),
crs = st_crs(trees_cd7)
)
# 3) Keep only density points inside District 7
cd7_poly <- cds |> filter(CounDist == 7)
inside <- st_within(density_pts, cd7_poly, sparse = FALSE)[, 1]
density_cd7 <- density_pts[inside, ]
# Pull back lon/lat for ggplot
density_cd7_df <- density_cd7 |>
mutate(
lon = st_coordinates(geometry)[, 1],
lat = st_coordinates(geometry)[, 2]
) |>
st_drop_geometry()
# Normalize density 0–1
density_cd7_df$density_norm <- scales::rescale(density_cd7_df$density)
# 4) Zoom box around District 7
bbox <- st_bbox(cd7_poly)
x_range <- bbox$xmax - bbox$xmin
y_range <- bbox$ymax - bbox$ymin # <-- this was the bug
x_lim <- c(bbox$xmin - 0.1 * x_range,
bbox$xmax + 0.1 * x_range)
y_lim <- c(bbox$ymin - 0.1 * y_range,
bbox$ymax + 0.1 * y_range)
# 5) Plot
ggplot() +
# NYC council districts background
geom_sf(data = cds, fill = "grey95", color = "grey80") +
# density tiles only where points are inside CD 7
geom_raster(
data = density_cd7_df,
aes(x = lon, y = lat, fill = density_norm),
interpolate = TRUE
) +
# outline for District 7
geom_sf(
data = cd7_poly,
fill = NA,
color = "black",
linewidth = 1
) +
scale_fill_gradientn(
colours = c(
"#2d004b",
"#54278f",
"#756bb1",
"#f768a1",
"#fe9929",
"#ffff99"
),
name = "Tree density\n(relative)"
) +
coord_sf(xlim = x_lim, ylim = y_lim, expand = FALSE) +
labs(
title = "Tree Density in NYC Council District 7",
subtitle = "Zoomed-in heatmap, masked to district boundary",
x = "Longitude",
y = "Latitude"
) +
theme_minimal() +
theme(
panel.background = element_rect(fill = "white", color = NA),
plot.background = element_rect(fill = "white", color = NA),
legend.position = "right"
)